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SUMMARY 

The general method for analyzing steady subsonic potential 
aerodynamic flow around a lifting body having arbitrary shape 
is presented. By using the Green function method, an integral 
representation for the potential is obtained. Under small 
perturbation assumption, the potential at any point, P, in the 
field depends only upon the values of the potential and its 
normal derivative on the surface, Z e , of the body. Hence if 
the point P approaches the surface of the body, the represen- 
tation reduces to an integral equation relating the potential 
and its normal derivative (which is known from the boundary 
conditions) on the surface (two-dimensional Fredholm integral 
equation of second-type ) . The question of uniqueness is exa- 
mined and it is shown that, for thin wings, the operator be- 
comes singular as the thickness approaches zero. This fact may 
yield numerical problems for very thin wings. However, numeri- 
cal results obtained for a rectangular wing in subsonic flow 
show that these problems do not appear even for thickness ratio 
t = .001. Comparison with existing results show that the 
proposed method is at least as fast and accurate as the lifting 
surface theories. 
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SECTION 1 


INTRODUCTION 

A general method for compressible potential aero- 
dynamics has been developed by one of the present authors*. 
Presented here are applications of the method to linearized 
steady subsonic flow. The method and the associated 
numerical procedure for evaluating surface pressure are 
described for general complex aircraft configurations. 
Numerical examples for finite-thickness wings are included. 
1.1 Definition of Problem 

The importance of an accurate evaluation of the flow 
field and related pressure distribution around an aircraft 
is well known to the aerospace designer. In the methods 
usually employed for the evaluation of the pressure 
distribution over a wing, the wing is assumed to have 
zero-thickness (.lifting surface theories, see for instance 
Refs. 1 to 4) . More recently more sophisticated techniques 
have been introduced for the evaluation of wing-body 
combinations. An excellent review of the state of the art 
in this field is given by Ashley and Rodden (Ref. 5) and 
therefore is not repeated here. It may be noted however 
that, as pointed out in Ref. 5, these newly developed 
techniques "are really just another evolutionary step 
toward the remote goal of flow-field analysis for wholly 
arbitrary configurations". On the other hand, arbitrary 


* See Ref. 6 




configurations have to be considered for aircraft design. 
Therefore methods for completely arbitrary configurations 
are becoming more and more important. A general method, 
suitable for use in automated design is presented here. 

A brief summary of this method is presented in the rest 
of this Subsection. The mathematical formulation of the 
problem is outlined in Subsection 1.2, while an outline 
of the report is given in Subsection 1.3. 

A general theory of the unsteady compressible potential 
flow around lifting bodies having arbitrary configurations 
and motions is presented in Ref. 6. In particular, for 
steady subsonic flow, the theory yields an integral 
equation relating the potential on the surfaces of the 
body to the "normal-wash”. However, the integral equation 
becomes singular as wing thickness goes to zero. Hence, in 
order to demonstrate the feasibility of the method, the 
problem of a rectangular thin wing is solved numerically 
herein. The method used consists of dividing the surface 
of the wing into N small elements, and approximating 
each element with its tangent plane at the centroid of 
the element (discrete element method). Then, assuming 
CJ> to be constant within each element and satisfying the 
integral equation at the centroid of the element, yields 
a system of N linear equations with N unknowns, the 
values, cp. , of the potential at the centroids of the 
elements. The results indicate that the method can be 
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used even for very small thickness (i.e., thickness ratio 
T =.1%. Loss of significant figures becomes important only 
when x decreases below about .01%. However, even for 
T =,01% the final results are in agreement (to plotting 
accuracy) with the existing ones. 

It may be noted, however, that for thick wings and bodies, 
the tangent-plane approximation may not be an accurate one. 
Furthermore, the discrete element method is not necessarily 
the best method. These questions are analyzed in Ref. 7 
where different numerical procedures are considered. These 
include a "modal method" (as an alternative to the discrete 
element method) whereby the unknown function tj> is expressed 
as a linear combination of prescribed functions (modes) with 
unknown coefficients. The results indicate that the discrete 
element method is preferable for both accuracy and computer 
time. Variations of the discrete element method are also 
considered in Ref. 7. These include a purely numerical co- 
efficient (without tangent-plane approximation) as well as a 
correction (.through numerical evaluation of the difference) 
of the error introduced with the tangent-plane approximation. 
Results indicate that the last procedure is the most accurate 
one. However, the improved accuracy does not compensate for 
the increase in computer time. In other words, if the com- 
parison is made with equal computer time (rather than with 
equal number of unknowns) the pure tangent-plane approximation 
is more accurate (for T =10%) than the corrected tangent- 
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plane approximation. Hence, the former one is used here 
although the second one is also mentioned for the sake of 
completeness. 


1.2 Formulation of the Problem 

In the following, the isentropic inviscid flow of a 
perfect gas, initially irrotational , is considered. Under 
this hypothesis, the flow can be described by the velocity 
potential $ . Consider a frame of reference such that the 
undisturbed flow has velocity in the direction of the 

positive x-axis. Then, it is convenient to introduce the 
perturbation potential , such that 

(pc+f) (1 ’ 1) 

Note that 2 0 in the undisturbed flow. The equation for 
the linearized potential subsonic aerodynamic flow is, 



( 1 . 2 ) 


where V 2 is the Laplacian in x,y,z variables, and 

M = (1 - 

is the Mach number. 

A very general approach is considered here by assuming 
that the body immersed in this flow has arbitrary shape. 
Thus, the surface of the body is represented in the general 
form. 


where the subscript B stands for Body. The boundary con- 


dition on the body is given by, 

-0 


(1 


By using Eq. (1.1), Eq. (1.5) yields. 


v f ■ v5 b = - 3 Se/ax 


or 


with 


W = LP 


^ " |V5J 3 


35, 


B 




(1 


(1 


(1 


Furthermore, as mentioned above, the boundary condition at 
infinity is given by 

Lj> 2 0 (1 

Finally, the pressure can be evaluated from the Bernoulli 
Theorem, 




(1 


or the linearized Bernoulli Theorem, 


, 2 a^p 


= - P n' 

'®0 'oc W «» 


which yields, for the pressure coefficient. 


c., = 


t-v„ 






3% 


(1 


(1 


.5) 

. 6 ) 

.7) 

. 8 ) 

.9) 

. 10 ) 

. 11 ) 

. 12 ) 


5 


1.3 ^Outline of "Report 


The method of solution presented here is based upon 
the well known Green function technique. The Green function 
for the subsonic flow is used in Section 2 to derive an 
integral representation of the potential cj> at any point 
in the field (control point) in terms of the values of 
and /dn on a surface surrounding the body and the wake. 

It is shown in Section 2 that the integral representation 
can be simplified under the assumption of small perturbation. 

In Section 3, the problem of a finite thickness wing 
in subsonic flow is analyzed within the limits of small- 
perturbation assumption, that is, linearized potential 
equation, boundary condition, and Bernoulli's theorem 
are used. Note that the surface of the body is not 
modified in view of the fact that the method is intended 
for complex configuration analysis. It is shown (Appendix A) 
that if the control point is on the surface of the body, the 
problem reduces to an integral equation. The question of 
existence and uniqueness is also discussed in Section 3. In 
Section 4, the numerical procedure for solving the integral 
equation is presented. In particular, it is shown that, 
for a zero-thickness flat wing, the operator becomes singular 

Hence, in order to verify the limit of applicability for 
low values of the thickness ratio, the subsonic flow around 
thin wings is solved numerically. The results are presented 
in Section 5, where the conclusions are discussed also. 


- 6 - 


SECTION 2 


INTEGRAL REPRESENTATION FOR THE VELOCITY POTENTIAL 

IN SUBSONIC FLOW 


2.1 Prandtl-Glauert Transformation 

In order to obtain a formally simpler expression for 
Eq . Cl .2), the following Prandtl-Glauert Transformation is 
made , 

*/p Y=# Z =$■ (2-D 

where 

<3 = ( \-fA z y 2 ( 2 . 2 ) 

Using these new coordinates, Eq. (1.2) can be rewritten as. 


r? 2 I r\ — Q 


V 

where V is defined as. 


V = 




sx 


(2.3) 


(2.4) 


Thus Eq. (1.2) is reduced to a Laplace equation. Eq. (2.3) 
is solved by the well known Green function technique, which 
will briefly be shown in the following subsection. 


2.2 Green's Theorem for the Laplace Equation 

The Green's function for Eq . (2.3) is defined as a 

function satisfying the following equation, 

= S (X-x-,, Y- Y, , Z-Z,) 12 . 5 ) 
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where S is the Dirac delta function. 

Multiplying Eq. (2.3) by G and subtracting from it Eq. 

(2.5) multiplied by ^ , one obtains, 

Gv 2 ?- <f v*<q = -fS ( 2 . 6 > 

Making use of the identity, 

v- ( a v b ) = ya-vb -f a v 2 b (2,7) 

one obtains, 

V‘( ( ^y ( f-fV^)=-fS (2.8) 

Defining a field function E such that. 


E 


and multiplying Eq. 
one obtains , 


1 in V (outside of body) 

0 otherwise 

(2.8) by E, then integrating it over V, 


(2.9) 


Iff 5 V; ( $ v, <f - y v; $ ) dv, = -Iff B f 8 d V, 

v V 

Using the following identities, 

Iff V' K, = <#> u . n ]dz, 

V J) 

Iff aS(x-x,,y-r,,2-2,;dv, = a ix, y,z) 

V 

with 



( 2 . 10 ) 


( 2 . 11 ) 

( 2 . 12 ) 


(2.13) 
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I 


o 

and noting that the solution to Eq. (2.5) is. 


4TJR 

with 

R = [(x-x, f + tY-Y.f+U-z ,) 2 \ 

Eq. (2.10) then becomes. 



(2.14) 


(2.15) 


4*E? = # 


where 



(2.16) 


(2.17) 


Eq. (2.16) is the desired integral representation for the 
velocity potential . 


2.3 Small Perturbation Hypothesis 

In this subsection, considerable simplification of Eq. 
(2.16) is obtained by making use of the small perturbation 
hypothesis. It should be noted that the small perturbation 
hypothesis has been evoked by assuming that the governing 
equation of motion is represented by the linearized equation 
of the velocity potential as shown in Eq. (1.2). For neg- 
lecting of nonlinear terms implies that the derivatives of 
(p are much smaller than one, i.e., in compact form, 


-2- - 0 (e)« I 


(2.18) 


I 


9 


where b is a characteristic length of the problem. 

This hypothesis is now used to simplify Eq. (2.16) into, 


2 


where ^ is defined in Eq. (1.7). 

In order to do this, consider the boundary condition 
given by Eq. (1.6). Under the assumption of small perturbation. 


Vk r °<*)«\ 


or, according to Eq. (1.8), 


I 3S 

IvSI ?>% 


_ 0 (e)«\ 


( 2 . 20 ) 


( 2 . 21 ) 


This means that the x-component of the normal is small com- 
pared to the other two, hence 


IvSI = t VS I 


( 2 . 22 ) 


and 


(O _ yf'VS _ 7jf-V3 h 2 as I 
1 IVSI | VS | | VS | (2.23) 

S _7^VS 

IVSI 7- 

Hence, Eq. (2.16) can be approximated by Eq. (2.19). 
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SECTION 3 


NUMERICAL FORMULATION 


3 ♦ 1 Introduction 

In the preceding sections, a brief discussion of the 
theory of steady subsonic flow is given. In Section 1, the 
problem is formulated. In Section 2, the Green Function 
technique is applied to obtain the integral representation 
for the velocity potential , i.e., ^ at any field point 

is related to *-P and on the surface of the body and the 

wake. A simplified expression is then obtained under the 
hypothesis of small perturbations. 

It is obvious that the general formulation presented 
here has no closed-form solution except for a few very 
special cases. Hence, in general, the use of high speed 
computers will be required. Thus, the numerical solution of 
the problem as formulated in Subsection 2.3 (small pertur- 
bation flow around an arbitrary wing) is discussed here. 

It should be mentioned that the numerical problems arise 
especially on the treatment of a thin wing (see Subsection 
4.6). Thus, in the discussion, it will be assumed that the 
body under consideration is a thin wing, although the for- 
mulation is valid for any body with sharp trailing edges 
(.see Subsection 3.3). 
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3.2 Integral Equation Formulation; Existence and Uniqueness 
of the Solution 


The integral representation of the potential for steady 
compressible flow in Prandtl-Glauert variables is given by 
(see Eq. 2.19), 

*JiE<f = -#y n - ? S F ciz + #?sa ,<7f )dz: (3.i, 

X z 

In order to analyze the question of uniqueness of the solution, 
the surface X is replaced by a smooth surface X surrounding 
Cat very small, but finite, distance; the boundary layer thick- 
ness, for instance) the body and the wake (Fig. 2a). The 
wake is truncated at a very large, but finite, distance from 
the wing. If the control point is on the surface , the 

function E assumes the value 1/2 and Eq. (3.1) reduces to,* 

= (3 ‘ 2) 
Z' 2 

If the geometry of the wake is known and if is replaced 

by the value (see Eq. 1.8 ), 

% = H = - BS /3>X, / /|V,S| 13-3) 

which assumes on the body and the wake, then Eq. (3.2) is an 

* 

See Appendix A 
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integral equation relating the downwash integral to the 
unknown value of CjJ on the surface. Note that Eq. (3.2) 
gives the solution of the exterior Neumann problem and, in 
this case, the solution of the equation exists and is 
unique* > for any smooth (Lyapunov) surface (Ref. 9 , pp. 
620-621) . 

Next, the surface 2j (surrounding the body and the 
wake) is replaced by the surface 21 r composed of two branches, 
the surface 21 3 of the body and the surface Z*/ of the wake 
(Fig. 2b). Thus, Eq. (3.2), combined with Eq . (3.3), reduces 

to** 


2TJ f 




(3.4) 


+ // ( fu - % ) { R )dz 

„ W 

In Eq. (3.4), U is the upper side of the surface of the wake 
and hence, the normal is understood to be the upper normal. 


* 

Note that, for the interior Neumann problem, the solution of 
the equation is not unique, for any arbitrary constant can 
be added to the solution. Physically speaking, one might 
say that, for the exterior problem, this arbitrariness is 

eliminated by the condition y = 0 at infinity. 

** 

The source integral on the wake is equal to zero. 
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Note that only the value of on the surface of the body 

(not on the wake) is necessary to solve the integral Equation 
(3.4) . 

From physical considerations, the solution of Eq. (3.4) 
is "very close" to the one of Eq. 3.2. Thus, it will be 
assumed that, if the geometry of the wake is known, the 
solution of Eq. (3.4) exists and is unique. It should be 
emphasized however, that this conclusion is based upon phy- 
sical reasoning. However, this reasoning is questionable, 
as shown by the remarks given in Subsection 4.6. Hence, a 
rigorous mathematical proof of the existence and uniqueness 
of the solution of Eq. (3.4) would be highly desirable. 

However, there are still two important questions to be 
considered: first, the geometry of the wake and second, the 

special behaviour of Eq. (3.4) when the thickness of the wing 
goes to zero. These two questions are discussed in Subsection 
3.3 and 4.6, respectively. 

3 . 3 The Wake 

As mentioned in the preceding subsection, the surface 
of the wake in Eq. (3.4) is not known. Thus, Eq. (3.4), 
which is satisfied on the body, must be completed by the 
equation on the wake, which says that the velocity on the 
wake is tangent to the surface of the wake. Thus, one ob- 
tains two coupled integral equations, one on the body and one 
on the wake, with ^ unknown on the body and *f n unknown on 


14 



the wake , whereas is known on the body and 

is constant along the x-direction on the wake, since no 
pressure discontinuity can exist across the wake: 




(3.5) 


According to the Kutta condition, this constant value 
is equal to the value of at the trailing edge. Given the 

velocity on the wake, the geometry of the wake is obtained 
by the condition that the velocity is tangent to the wake. 
This approach has been successfully used in Ref. 10 to study 
the transient incompressible flow around a wing after a 
sudden start. However, from a practical point of view, this 
approach is too lengthy and a simplified treatment of the 
contribution of the wake is presented in the following. 

Note first that 


f&M )dz= -fM dz 


where 



(3.6) 


dz^ = dz c&d. (N* ) 


(3.7) 


is the projected area (along the direction ) and 

- dj-j" (3.8) 

Ft* 

is the solid angle (see Fig. 3) . 

Next, consider the wake integral as a sum of M strips 
in the x direction (see Fig. 4). Applying the mean value 
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theorem, one obtains (note that 4 ^ is only the function of 

Y, ) / 

A f^)^l-^)dz=-ff Af(X)dn 


u 


= -Z A f(Y,Jffdn = - z Afixjn 

fa-l Ti w m-l •* 


(3.9) 


rf\ 


M 

dn = - 

tC w=/ 

where A.(f>(Y irn ) are the mean values of 4^ for each strip 
Hi , and jQ m are the solid angles of the strip L\Yy\ • Eq. 
(3.9) shows that any changes of the wake such that solid 
angles -Q^are not altered, do not have any influence on the 
value of the wake integral I w . 

This suggests that a "reasonable" geometry for the wake 
can be assumed, provided that the values of the associated 
solid angles are not excessively different from the true ones. 
Hence, it is possible (and convenient) to approximate the wake 

by straight vortex-lines, parallel to the direction of the 

* 

flow, emanating from the trailing edge of the wing. For, 
geometrical considerations show that the solid angles, -O.^, 
are changed only slightly. With this assumption, the wake 
integral simplifies considerably and its contribution reduces 
to a line integral. For if the trailing edge is given by 

X =■ X t£. ( Y ) 

te.(Y) (3 ' 10) 

then the equation of the surface of the wake is given by 


5 = z, -z T . E .(t) = 


0 


(3.11) 


^This is consistent with the small-perturbation assumption. 

The limiting behaviour as the thickness of the wing approaches 
zero is discussed in Ref. 6, Appendix G (See also Subsection 
4,7 of this reports. 
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and 


b/2 


-<XS 


I.'-I C, Y, / Af<v,S-R)-fiidx l dr, 

' 2 X T E. 

= i / /V ,J wdX 

where b is the span of the wing and 


(3.12) 


J« = [(z,-Z)- ~'(YrY>] 


*TE. ^ 


- I 


(3.13) 


fl-( »] 

(y,-Y)+(Z r Z) zL R J1 


with Z, — In particular, if the trailing edge is 

in the plane Z t = 0 (i.e., Zj.JX)* 0) , Eq. (3.13) reduces to 



(YrY)-t-Z 




^TE. * 

[(x TE -/) 2 +(>;-rA(Zrz; 2 J 




} 


(3.14) 


In conclusion, .under the reasonable assumption of cylindrical 
wake (straight vortex-lines) the effect of the wake sim- 
plifies considerably and Eq . (3.4) reduces to, 


b/2 

* J " d 

with given by Eq . (3.13). 

Finally, an important remark about bodies • without sharp 
trailing edges, must be made. In the discussion presented in 


17 



this subsection, it was assumed that the wing had a sharp 
trailing edge. Note that the results can be easily genera- 
lized to the case of general bodies with sharp trailing edges. 
However, for bodies without sharp trailing edges, the inviscid 
flow theory is incapable, in general, of predicting the lo- 
cation of the stagnation point from which the wake emanates. 
This can be easily seen in the case of rotating cylinder of 
finite length. From the experiments the location of the 
stagnation point depends upon the angular velocity, cO . On 
the other hand, the equation of the geometry of the cylinder 
does not depend upon 60 and thus cd does not even appear in 
the equation of the inviscid flow. 

In the following, it is assumed that the body under 
consideration has a sharp trailing edge. However, a viscous 
theory which predicts the location of the stagnation point 
from which the wake emanates, is necessary in order to extend 
this method to bodies without sharp trailing edges. Similar 
consideration holds for the case of detached flow (when this 
can be approximated by a wake emanating from a point different 
from the sharp trailing edge) . 
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SECTION 4 


NUMERICAL PROCEDURE 

4.1 Discrete Element Formulation 

The governing equation for the velocity potential, after 
Prandtl-Glauert transformation, for steady subsonic flow, is 
shown in Eq. (2.19) and is rewritten here, 

4T e'f=-#%-kdz+#£j,(-*)'fdzi (4.i) 

r. X 

where ^ on the left hand side, represents the velocity 
potential at a field point P and ^ on the right hand side 
represents the potential at a point P^ on the boundaries 
(body surfaces and wake). E is defined as, 

E = 0 ifPis inside the body 

= 1/2 if P is on the surface (4.2) 

=1 if P is outside of the body 

From Eq. (4.1), it can be seen that if (j> on the boundaries 
is known, the velocity potential at any field point can be 
obtained immediately. So the first step is to find (j? on 
the boundaries. With E = 1/2, Eq. (4.1) becomes, 

2Jry = - # f*ik dz 

Ze 

~ , (4.3) 

+ //* 

TJ W ‘ 

_ T ,w 

where 2T B is the surface of the lifting body, U is the 

upper surface of the wake and. 
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(4.4) 


^ f = - t~r) 

on the wake. Eq. (4.3) cannot be solved analytically except 
for some special cases. Hence, numerical scheme has to be 
employed. The procedure used in this study is as follows. 

As mentioned in Section 3, the wake is approximated by a 
plane parallel to the direction of flow, and Ay on the 
wake is equal to Acf at the trailing edge and integration 

in the x-direction from x = x _ to x->«3 is carried out as 

t . b . 

outlined in Section 3.3. Secondly, the surface of the body 
is divided into discrete surface elements (see Fig. 4). In 
each element, ^ is approximated by its value at the centroid 
of the element, then the integration is carried out within 
each element. The resulting equation can be written as 


with 



-■Ai 


where 


= # 


( 


2T\Hj 


)dzc 



= [ ( x^-x, )\ ( Y 



(4.5) 


(4.6) 

(4.7) 


(4.8) 


is the distance between the dummy point of integration P to 

( k ) 

the center point P of the element k, and 
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for the boxes in contact with the trailing edge and 

Hi = 0 

otherwise. Eg. (4.5) can be rewritten as 


(4.10) 


N 


with 


ft 




(4.11) 


(4.12) 


where is the Kronecker delta. Solving Eq. (4.11), one 
obtains the solution of the problem. So the question which 
remains now is how to evaluate b^, c^ , and w^^ . 

For a wing symmetric about the plane y = 0, Eq. (4.10) 
can be rewritten as 


A 


A 

N 


2 


(.4.13) 


where N 1 = N/2 is the total number of elements on the right 
hand part of the wing and c^ is the influence of two elements 
(in symmetric position with respect to the plane Y = 0) on 
an element on the right hand part of the wing. 


Su 


5T fJA 


with 


YrV + If 1 ZrZ> ] 


(4.15a) 
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(4.15b) 


[ € (x - x; + 


where 



~ [ (x-x/i ty ,~Y) 2 -t <z,-zf ]' /z 

{L) 2 2 2 *-iV? 

- [(*,-*) -+ f y, + y; + iz,-z) ] 

Similar expression holds for w,. . 


Evaluation of c^, w^ and is shown as follows. 


4.2 The Geometry of the Wing 

The planform of the wing is defined as 

*n,= 

^.E, “ ^.E. ( 

i for o<ty < 0/2 


(4.17) 

(4.18) 


where b is the span of the wing. The chord of the wing at 
any point is then defined as 


C( ^ = 


(4.19) 


This planform can be transformed into a rectangular one by 
the following transformation 



X ~ %UE. 

X-TE. ~ %L.E. 


0 < £ £ 1 


(4.20) 


Z '- 2 V/b 



(4.21) 
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The thickness, h, of the wing is defined as 


4 = o-3)Ji-i‘ 


(4.22) 


where c is the maximum chord length and t is the thickness 
max 

ratio of the wing, i.e., 


t = 


Tr*sx. 




(4.23) 


Summing up, the geometry of the wing is represented by 


^ ^ ) J *+ 

5 = yi/z 

■for 0 £]£ | -/•<?$, 


(4.24) 


If the angle of attack cK, is different from zero, then 

X - X CCKL o( -t 5 

y = y (4.25) 

£ = CO<lo( 

For small £ and cK / Eq. (4.25) can be approximated by 

X - X 

}f = V (4.26) 

3-= I -%■=< 


As mentioned before, in the process of solving Eq. (4.4), 
the wing is divided into elements and the mean value of in- 
side each element is approximated by that at the centroid of 
the element. Because ^ and the slopes of the wing vary more 
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rapidly near the leading edge and the tips, smaller elements 
are required at these regions. Using the following trans- 
formation. 


5 = x 2 

*1- 1 - ( i-y) 2 


(4.27) 


then, dividing the transformed planform into even mesh (i.e., 
the mesh size is 


AX. = I /MX 
AY = I/NY 


(4.28) 


where N X and NY are the numbers of elements in the x and y 
directions, respectively), the boundary of the element on 
the n x row and n column is 


( n x -i) Ax 
( ny-i) AY & ri 


( 4 . 29 ) 


The centroid of this element is defined as 


x c = (n x --s) ax 

r t =(nj-^f (4 ' 

Note that Eqs. (4.28) and (4.29) yield smaller elements near 
the leading edge and tips . 


4.3 Calculation of c 


ki 


The evaluation of c^ consists of two parts. The first 
part is to approximate the wing by replacing each element by 
a plane tangent to the wing at the centroid of the element 
defined in Eq. (4.30). The boundary of the tangent plane is 
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defined by Eq. (4.29). The equation of the tangent plane can 
be written as 

3-h =(§i-) c (*-Zc)+(-5§)cW-Vc) 

where the subscript "c" indicates the value at the centroid 
mentioned before. 

The integration of Eq. (4.14) on these tangent planes 
can be carried out analytically as shown in Appendix B. 

This tangent plane approximation yields good results for 
thin wing (thickness ratio varying from .1% to l.%). The 
computer time required for this tangent plane approximation 
ranges from 13 sec. for NX = NY = 4 to 128 sec. for NX = 

NY = 7 , on the IBM 360/50 of the Boston University computing 
center. 

For a thick wing, the tangent plane approximation may 
not yield good results. Then the solution can be improved in 
the following manner. The most severe error comes from two 
sources. First, the curvature of the surface changes very 
rapidly near the leading edge and the tips. So, the tangent 
plane is not a good representation of the true surface. Se- 
condly, as can be seen from Eq. (4.16), ^ 0 when k = 

i (i.e., field point coincides with source point), so 1/R^ 
in Eq. (4.14) varies rapidly, hence again making the tangent 
plane a poor approximation. Thus, it is necessary to inte- 
grate numerically the difference between the actual surface 
of the element and the tangent plane; and add this correction 
value to that of the tangent plane approximation. 
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The numerical integration scheme used here is the Gaussian 
quadrature. For integration over the leading edge and the 
tips (but i# j), a four point Gaussian quadrature is used. 

For i = j, a polar Guassian quadrature (see Appendix C) is 
used. The results obtained in Ref. 7 indicate that no correction 
is needed for the elements which are not on the leading edge nor 
on tips and are such that i j . 


4.4 Calculation of w^ 

As mentioned before, the surface of the wake is approxi- 
mated by a surface parallel to the direction of the flow, 
emanating from the trailing edge of the wing. Based on this 
assumption, it is shown that Eq. (3.12), 


with 


'2 


(4.32) 


(4.33) 


■ ri2, f 2 ,-^ t-m i- J 


and ^ • Then, integrate Eq. (4.32) and 

lump to the last row of elements in contact with the 

trailing edge. It should be noted that if dZ/dY=0, or if 
this term is negligible compared to the other term in the 
same bracket, Eq. (4.32) can be integrated exactly. 
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4.5 Calculation of 

The process of evaluating b k , i.e., Eq. (4.6), is similar 
to that for evaluating c k ^ . For thin and moderately thin wing, 
b^ is evaluated analytically on the tangent plane approximation. 
Integration for b^ on these planes is shown in Appendix B. For 
a thick wing, the solution is improved, as for c^ , by numeri- 
cally integrating the difference between the real surface and . 
the tangent plane and adding this difference to the value ob- 
tained by the tangent plane approximation. 

It is worth noting that, for a wing with a plane mid- 
surface (symmetric about this plane) , |f n consists of the 
symmetric part (value on upper surface is same as that on 


lower surface) and the antisymmetric part. The symmetric part 

fn S> 13 ' 


I 

in 1 ** 

and the antisymmetric part is, 

W 

% 


(4.34) 


(4.35) 


where the upper sign stands for upper surface and the lower 
sign for lower surface. Since b^ is linearly proportional to 

, b fc can be separated into symmetric part and antisymmetric 
part. Eq. (4.11) can be rewritten in two equations, 



(4.36) 

(4.37) 
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( S ) (*) 

where ji stands for the symnietric solution and for anti- 
symmetric solution. Separation of into these two parts 

avoids elimination of significant figure due to the presence 
of the other, hence greatly improves the accuracy of the so- 
lution. 


4.6 Calculation of the Pressure Coefficient C 

E 

As shown in Eq. (1.12), the pressure coefficient can be 
written as 

. “2 d< f/ax (.4.38) 

As a first approximation, C^ is evaluated by finite difference 
method , i .e . , 


( c p ) i - -z ( fj+r ft ) _ _ 2 ViH-fi 


(.4.39) 


AX; ^ X 2 X, c 

where (C ) . is the C at the rear end of the element on the 

pip 

ith row, ft is the value of ^ at the centroid of the element 
on the ith row and is the distance, in the x-direction, 

between the centroids of the ith and i + 1st rows of elements. 

A better approximation for C^ is to combine Eqs. (4.1) 
and (4.38) and calculate C^ by integration. For the present 
purpose, however, Eq. (4.39) is a good approximation. 


4.7 Limiting Behavior for Zero Thickness 

As mentioned above, the formulation described thus far 
becomes singular in the case of zero thickness. This is 


28 


shown clearly by the fact that, for lifting surface (in the 
plane = 0) , Eq. (4.1) reduces to 


and 

= ffj K-%) 4 f -k) zr J x ' d Y ' 


14.40) 


C4.41) 


where Z, is the portion of the plane = 0 (upper side) which 
contains the wing and the wake. By adding and subtracting Eqs. 
C4.4 q) and (4.41), one obtains 


and 


C4.42). 



C4.43) 


This implies that (since, as well known, there exists a non- 
trivial solution the operator shown in Eq. (4.43) is 

singular . 

Hence, one can expect that the numerical procedure also 

has a singular behavior. In order to show that this is indeed 

the case, consider a symmetric wing with angle of attack oC 
and thickness ratio T , and let go to zero. In this case, 

Eq. (4.6) shows that 

i £x> 1 ’a = 0 (4 - 44) 


since (^) u In order to simplify the discussion. 



the numbering of the elements is assumed to be such that the 
odd (even) numbers correspond to elements in the upper (lower) 
surface and that the element in opposite position to the 
upper element i, has the number i + 1 (see Fig. 4). As 
before, upper element i and lower element i + 1 will be 
called "opposite elements". 

With this numbering, it is easy to show 6 that, according 
to Eg. (4.7) , 


liv* [Cfc] = 


o 


-f I O O I 
1 I 

O I o O |_ 


0 O 1 O -I I 

° 0 JLT 1 0 ± 
i I 


[ o o 
\0 o 



o 

CP 


i 


(4.45) 


O O 1 0 o 
I 

0 0 \ o O 


I o -/ 


In other words, all the coefficients c^ are equal to zero 
except for the ones relating opposite elements, which assume 
the value -1. Furthermore, the coefficients w^ are equal 
to zero. Hence, Eq. (4.11) in the limit, as £ goes to zero, 
reduces to 


I I I o o i 

_o_[_ 

o o i I f i 

o o l i (I 


L 


O o j o o j 

O o I O o I 

i i 


o o 

o O 


o o 
O O 


H 1 

!' - 

! i I 


? 4 =0 


(4.46) 
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This equation can have a nontrivial solution since the deter- 
minant is equal to zero. 

Note that this result implies that zero thickness wings 
(lifting surface theory) are more difficult to deal with 
than finite thickness wings. 

However, this shows also that, by using the method pro- 
posed here, one may encounter numerical complication due to 
the fact that, for very thin wings, the determinant is close 
to zero and hence, one may encounter strong elimination of 
significant figures. This implies that one has to be very 
careful in the evaluation of the coefficients c^ and b^. 

In order to establish the practical limits of the 
applicability of the method, Eq. (4.11) has been solved 
numerically for very small values of £ . The results are 

presented in Section 5. 
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SECTION 5 


NUMERICAL RESULTS 


5.1 Introduction 

Some numerical results showing the validity and flexibility 
of the present method are shown in this section. The con- 
vergence of the present method is analyzed in Subsection 5.2.1. 
The effect of the numerical correction for a thick wing (as 
discussed in Subsections 4.2 and 4.4) is discussed in Subsection 

5.2.2. The effect of thickness is described in Subsection 

5.2.3. Finally, comparison with existing results is given 
in Subsection 5.3. 

To give an idea about the functional type of the solution, 
the value of ^ and C^ are shown in three-dimensional forms 
in Figs. 5 and 6, where 

c r' 4C rVV (5 - 1) 

is the lift coefficient. The thickness ratio is 2 : = .001, 
angle of attack cK. = 5°, aspect ratio b/c = 3 and the mesh 
NX = NY = 7 (that is N = 98) . It may be noted that the diagram 
of is flat in the neighborhood of the root (more precisely 
0 at root). Hence, in the following discussion 
of convergence, the value of at the centroid of the elements 
in contact with the root (root elements) is used. 


5.2 General Characteristics of the Method 

As mentioned above, in this subsection the convergence. 
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the effects of the numerical correction and of the thickness 
are analyzed. 

5.2.1 Convergence 

In order to study the convergence of the solutions, the 
case C = .001 was solved for NX = NY = 4 , 5, 6 and 7, res- 
pectively. The values of the lift distribution ^ C 

XT 

at the root elements are shown in Fig. 7. 

The results show that the solution is convergent very 
fast and that the case 4 x 4 is sufficient for an accurate 
analysis. The computer time employed on the IBM 360/50, 
available at the Boston University Computing Center, is 
given in Table 5.1.* 


Number 

of Elements 

TABLE 5.1 

Computing Time (Sec.) 

4 x 

4 x 

2 


13 

5 x 

5 x 

2 


30 

6 x 

6 x 

2 


64 

7 x 

7 x 

2 


128 


5.2.2 Effect of Numerical Correction on Thick Wing 
Table 5.2 shows the comparison between the results of 
the tangent plane approximation and that with numerical 
correction for a wing with thickness ratio £ = .2. The other 
data are the same as those for Fig. 5. The computer time 

* 

Advantage of symmetry with respect to Z was taken. 
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required for both solutions are also shown. It is seen that 
even for this thick wing, the tangent plane approximation 
yields good results and considerable computer time is saved. 
So all the discussion given in the following subsections is 
based on results from this tangent plane approximation. 


TABLE 5.2 



T.P.A. 

T.P.A. + N.I. 

C JL 

8.09 

7.91 

at 

2.60 

2.59 

root 

.975 

.900 

Computer 

Time 

(sec) 

22 

136 

5.2 

.3 Thickness Effect 



In order to analyze the thickness effect, the problem 
has been solved for four values of the thickness ratio, 

= .1, .01, .001 and .0001, respectively. In all these 

cases, the number of elements in both x and y directions is 
NX = NY = 4 . Hence, the total number of elements (for upper 
and lower side of the right half of the wing) is N = 32 
(i.e., Eq. 4.13 is a system of 32 equations and 32 unknowns). 
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For the value £ = .001, no message indicating strong 
elimination of figures was given, whereas, for the value 
£• = .0001, a message indicating an elimination of signi- 
ficant figures higher than the prescribed tolerance at the 
19th step was obtained.* 

Hence, only the cases £ = .1, .01, and .001 are pre- 
sented here. The lift distribution Co = - AC at the root 

-< p 

elements, is shown in Fig. 8. The results indicate that the 
solution converges to a zero-thickness solution and that the 
solution for = .001 is a good approximation for the zero 
thickness solution. It may be noted that the results here 
showed that a thinner wing yields a higher lift. 

Note that, according to Eq. (3.6), 


M 

Z 

C-) 


5 = 


C AL ' 


27f n 

Zl 




27 f 


(5.2) 


since the point from which the solid angle is evaluated is 
on the surface 2. . This equation is poorly satisfied on 

the leading edge and the tip where the approximation of 
the surface element with its tangent plane is poorer. The 


* 

For the solution of Eq. (4.13), the standard IBM SUBROUTINE 
GELG has been used. The value of the tolerance (which is 
compared to the ratio between the pivot at the nth step and 
the initial step) was chosen to be TOL = .001. 
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r 


N 

poorest values of 2 C. . at tip or leading-edge elements" 

C-l 

are given in Table 5.3, column 1, whereas the poorer value 
for the "internal element" (not at the leading edge nor at 
the tip) are given in column 2. 


TABLE 5.3 


X 

1 

2 

.1 

.76946 

.98986 

.01 

.96247 

.99899 

.001 

.99625 

.99990 

.0001 

.99960 

.99999 


Table 5.3 indicates that, for the case ^ = .1, the approxi- 
mation of the surface elements with its tangent plane is not 
satisfactory for the "non-internal elements". 

5.3 Comparison with Existing Results 

In order to evaluate the accuracy, comparisons with 
existing results are given in Figs. 9 to 14. Different 
values of Mach number, aspect ratio and different planforms 
are tested to show the flexibility of the program. Thickness 
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ratio shown in the figures are chosen arbitrarily, 
are defined as 


£ = X/Cft) 

\ -- 2 y /t 


^ and ^ 

(5.3) 


Figure 9 shows the comparison of the value between the 

11 12 

present solution and Cunningham's Theory and experiment 
on three cross sections of a rectangular wing. As can be seen, 
agreement between the two theories is very good and agreement 
between theory and experiment is very good except when near the 
leading edge. Figure 10 compares the results obtained with the 
present method and those obtained by Hsu (Ref. 13) , Kulakowski 
and Haskell (Ref. 14) and Cunningham (Ref. 11), for a rectangular 
wing of aspect ratio AK = 1, thickness ratio T= .001 and Mach 
number M = .2. It may be noted the excellent agreement between 
the present method and the ones of Refs. 11 and 14, which are 
generally considered to be very accurate. 

The results for a delta wing ( AZ =2.5, T = .005 and 
M = 0) are compared in Figure 11 with the Widnall results (see 
Fig. 7.7 in Ref. 15). The results for a tapered swept wing 
( /R. = 3.0, taper ratio = .5, A _l = 45°, M= 0.8, and oi = 
4.13°) are compared in Fig. 12 with Cunningham's results (Ref. 

11 ) . 


Finally, it should be noted that the comparison of the 
integrated values of the pressure is generally more signifi- 
cant than the comparison of the pressure distribution because 
small (relative) errors on the pressure distribution near 
the leading edge can yield large errors on the integrated 
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values. Hence, the values of the section lift coefficients 

/ 

1 Ac f« < 5 - 4 > 

o 

are compared in Figs. 13, 14 and 15 to the ones obtained by 
Yates (Ref. 16) for a rectangular wing with aspect ratio 
& = 4,4,8, thickness ratio t = .001 and Mach number 
M = 0, .507, 0, respectively. It may be noted that the 

Figs. 13 and 14 (incompressible and compressible flow, res- 
pectively) show an excellent agreement, whereas in Fig. 15, 
the agreement is less satisfactory. 

In conclusion, the agreement is excellent for the cases 
presented in Figs. 9, 10, 13, and 14, but less satisfactory 
in Figs. 11, 12 and 15. The reasons for the small disagree- 
ments shown in Figs. 11, 12 and 15 are now being analyzed. 
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SECTION 6 


CONCLUDING REMARKS 

A general method for analyzing steady subsonic flow 
around lifting bodies of arbitrary configuration is presented 
in the preceeding Sections. The method is then applied to thin 
wings in steady compressible flow, which is the most challenging 
problem for this method. The results indicate that the method 
is accurate and does not require excessive computer time. The 
rate of convergence is surprisingly high. Similar results 
were obtained for oscillating wings in subsonic flow (Ref. 18). 

Since the problem of elimination of significant figures 
is not encountered even for very thin wings (thickness ratio 
T = .1%) it is expected that no such problems will be 
encountered in applying the method to realistic complex 
configurations. Extensions of the method to complex configura- 
tions are now under consideration. 
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Fig. 1 Geometry of the problem 
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Fig. 5 Potential distribution for a rectangular wing with aspect ratio b/c = 3, 
angle of attack 5° and Mach number M = .24 
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Fig. 6 


Lifting pressure distribution for a rectangular wing with aspect 
ratio b/c = 3, angle of attack o(= 5° and Mach number M = .24 
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Lifting pressure coefficient at root boxes for thickness ratio 
T = ,001 and different values of NX and NY for a rectangular 
wing with aspect ratio b/c = 3, angle of attack °C= 5° and Mach 
number M = .24 




Effect of decreasing thickness on the numerical scheme: Lifting 

pressure coefficient at the root boxes for NX = NY = 4 and 
different values of the thickness ratio x , for a rectangular wing 
with aspect ratio b/c = 3, angle of attack oC = 5° and Mach number 
M = .24 





Lift distribution for a rectangular wing 
with aspect ratio b/c = 3 , angle of attack 
oi = 5° and Mach number M = .24 


Fig. 9 
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Fig. 10 Lift distribution per unit angle of attack 
at 2y/b = .7 for a rectangular wing with 
aspect ratio A.R. = 1 and Mach number M = .2 
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Fig, 12 Lift Distribution at 2y/b = .707 for a swept 
tapered wing with aspect ratio A.R. = 3, 
taper ratio T.R. = .5, A. , /4 .= 45°, angle 
of attack oC= 4.13® and Mach number M = .8 
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Fig. 15 Section liftrcurve slope for rectangular 
wing with aspect ratio A.R. = 8 and Mach 
number M = 0 




APPENDIX A 


THE VALUE OF THE FUNCTION E ON THE SURFACE 
A . 1 Introduction 

Equation (2.19) gives a representation of the potential 
anywhere in the volume, in terms of the values of and ^ 

on the surface of the body. The values of 'd^/'dfL on the 

surface of the body are given by the boundary conditions, but 

the values of are not known. In order to solve the pro- 
blem, it is thus necessary to obtain first the values of ^ 
on the surface. This can be done by letting the point P of 
the volume V approach a point P* of the surface. 

In this Appendix, it is shown that, in the limit, Eq. 

(2. 19) is still valid if the definition of the function E is 
generalized as follows 

E = 1 outside 21 

E = 1/2 on Z (A.l) 

E = 0 inside 21 

By letting P — > P* , the integrands become singular in the 
neighborhood of P* . Thus, it is convenient to separate the 
contribution of a small neighborhood* of P* , which will be 
indicated as Zg. . 

For steady compressible flow, Eq. (2.19) can be rewritten 


The neighborhood Zg. is a small circular surface element, 
with center P* and radius £ . 


58 



as 




z-z< 


z-r € 




(A. 2) 


where is the contribution of the neighborhood of P* , given 


by 


The analysis of &e is highly simplified by using local 
A A Z' A 

coordinate X, Y, Z, with Z normal to the tangent plane 
(directed from E = 0 to E = 1) . Then, separating terms of 
order £ , Eq. (A. 3) reduces to* 


V-^SUI 


I 


A 


Bz ' *x,*+Y , 2 <g 2 ^ X, 2 t f,*t ^ 3 


-% a 


dx, dr. 


c/A.dyi -j. o ce; 


(A. 4) 


VFty.Vz 2 

where the subscript * indicates evaluation at P* . By using 
polar coordinate 

R = 7 x, 2 + r, a 

e =■ An" x/x, 

one obtains 


(A. 5) 


* x\ ✓a -A 

Note that X = Y = 0, Z 1 


0 and 



Z,=0 
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^ = - z7 i ci 1 1 


' dZ , '* 


J 


A2, ^ 2 

RfZ 


S\ . y\ 

Rcl R. 


-*n 


J 


& 


ft * c ) z. r ^2 zs 2 

0 7 R t z -2 


/ A , 2 

1 RdF{ -t 0 (e) 


Noting that 


3 ( 


2 SZ "/^F 

Eq. (A. 6) becomes 




r a « y^i 5 


a<p 


& 6 - ~ 2y r ( J* [JW+P ] 0 

~ 27{ % * f/^zp ] a -i-O(e) 

y\ 

Finally, by letting P go to P* , (that is, Z— ► 0) , one 


l'" S e = iL" [ ' 27r ( §?)* ( I^ 2 -12 1) 


P-P A 


z-»o 


( 


7e 2 -f z 2 J z/ 


)]+ 0(e) 


-f] + o(e) 


A 


3Z, /* 


121 


- -t 


A^T?* + 0 (e) 


(A. 6) 


(A. 7) 


(A. 8) 


obtains 


(A. 9) 
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where the upper (lower) sign holds for Z>0 (Z<0), that is, 
when P originates outside (inside) the surface 21 / corres- 
pondingly, the function E assumes the value E = 1 (E = 0) . 
Finally, using this result in Eq. (A. 4), one obtains 


47J U=±i)% = -ff g+dZ+ff fir^dz + OCe) 

z-z e 


(A. 10) 




Note that, in both cases (P inside or outside 21 ) , 


E* = E ? Z = T = { 

= o ■*■£ - 2 

Furthermore, R* is the distance between the dummy point, 
P^, and the control point (on the surface 2. )» P* . Hence, by 
letting 6 go to zero, Eq . (A. 10) yields 


(a. 12 ) 


It should be emphasized that the limit 6"^ 0 is now per- 
formed with P on the surface 21 • This implies that the con- 
tribution of Zg is now of order £■ . In order to clarify 

this point, consider the quantity 



(A. 13) 
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and note that 


(im { lim { 


A 

2 


•A 

z 


Z ~>0 
whereas 


Z-5>0 


’Je 7 +z 2 


izl 




(A. 14) 


t f m { Jim ] =r I i’m 0=0 


e o -^0 


^0 


(A. 15) 


The difference between these two limits is due to the fact 


that, in the limit (as Z— >0), the integrand of I^behaves 
like a Dirac delta function and hence, its contribution for 
a domain which excludes the singular point is zero. 

It may be worth noting that, in Eq. (A. 9), the sequence 
of limits indicated in Eq. (A. 14) must be performed, whereas 
in Eq. (A. 12), the one indicated in Eq. (A. 15) must be used. 

Finally, it is shown that the results obtained here are 
equivalent to the definition of E given by Eq. (A.l). 

Note that Eq. (2.19) must be used if P is outside or 
inside the surface, whereas Eq. (A. 12) must be used if P is 
on the surface. However, by comparing Eqs. (A. 12) and (2.19), 
it is easily seen that Eq. (2.19) is valid everywhere (outside, 
inside and on the surface Z ) , if the convention is made 
that E is given by Eq. (A.l). 
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APPENDIX B 


SOURCES AND DOUBLETS ON A TRAPEZOIDAL ELEMENT 


B.l Trapezoidal Element 

In this Appendix, the effect of the sources and doublets 
distributed on a trapezoidal planar element, are obtained in 
analytic form. As mentioned in Section 4, it is of interest 
to consider planar elements described by the equation 

z,-z, c = * (x,-x, c ;tj3(Y,-Y, c ; (E - 1 

where the subscript c stands for centroid of the element, 
defined in Eq. (4.30). The boundary of the projection of 
this element on the plane Z = o is given by 

X lm ( Yi*" Y) ^ X| X\ip + ^ X) (B. 2 

X m * T, 6 Y,r 

Equation (B.2) represents a trapezoid and the element defined 
by Eqs. (B.l) and (B.2) is called trapezoidal planar element 
(see Fig . B.l). 

B.2 Doublet Distribution 

Consider first the integral of a doublet distribution 
of unit density over a trapezoidal planar element, given 
by 
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where 


ic> - II Zj, ( Ft ) dz (B.3) 

Xp Kp+dftYrW i 

= i^{ , d (r r , &kl ff 1 1 YrYh M (ZrZ)]dx ' 

A |rn ru m v iris 

r, T Vttyt-T), 

= S J dr.f L°< (X.-X)-t p (Y,-Y)~ (z rZj]-^dx, 

Y,m X, m td M (YrY) 


s = - 3S/az ' = , 


P s /aZ, | 


= -1 


on upper surface 
on lower surface 


(B. 4 ) 


and use has been made of the fact that, according to Eq. (B.l), 


Z = /*s 

a *i 7 az } 

pS/* 3 


(B. 5 ) 


* Y t / 

It should be noted that, for a trapezoidal planar element, 
^ and /3 are constants and that, according to Eq. (B.l), 

1= <2 r 2)-3.<x,-x) -p(Y,-Y) 

= (z,-z)-x(x icr x)-p ( YrY) 

Thus, Eq. (B.3) reduces to 


(B . 6 ) 


^ _ Ti-p Xi-ptdytYrY) 

l c -ssJdYj ~Z7 dx i 

Ylrn x m - td ^ YrY ) n 


(B.7) 
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with 


Vz 

r = {(x,-x) 2 -t ( Y,~yA [(§-+»< ip(X~Y)T} 

= £[(■+* 2 ) /2 ix,-xh*( p < k-y)+ §)]*+ (Hefrp 2 ) (B> 8 } 

t [ ( i+&p a /*(YrY)+pV( i +<* V ^f/a+o?)} A 2 

= f j 3 ] /2 

where 

3 =( i+^/l x-x) -M [ p(Y~Y) + s ] / 0+^ J ) 

£ = [(i+S^^'^Y.-Y) •+ p&Ah&p^T/i iiZ ! ) ' /2 (B - 91 

f = § 


Integrating Eq. (B.7), one obtains 


Yr 


f _ -ss f r r_i±_ I 

D f +SW. 1 «a. ,J2. 2 2C2,!4J“" 


I, 




l d l 


(B. 10) 


where 


S f = f / 1 * 2 / 2 f X,^X*d p ( VY>] +* t P ( Y~ YHS ]/[/ + a 1 ] 


"2 

Cb. ll) 


Jo-b ** ^ 


and 


|^= ( t p (y / ,-y)+ s]/t'-K* 2 ) ,/2 


CB.12) 



with 


Xf ( 1 + * Z ) [ x if -X + 0* - ol p p ) S' /( I + *2 + jg*)J 

lif- [off + (1 + ^ 2 )^ P ] /ft + +p 2 j ,/2 

(B . 13 ) 

£* ■ - ( 1 1 • : 1 )" 2 [x. *, - V + (* - 4 * £ ) 5 /(! +**+?)] 


. I,„ - ^4 (1 4 W^]/[/ + Z*+p z ] ,/2 

a 

Substituting Eqs . (B.ll) and (B.12) into (B.10) for and 

and changing the variable of integration from Y, to * 
one obtains 


* 


T - -5S 

-*-h , . —2 — 


A 

.V 


Jr 


7«* Y Jj 


rp 


A 

7 


3om"t l™ 1 ? 


(l+52+ ^Y ( ? + 5 4 ;r * 


ft 


(B. 14 ) 


where 


? w = Em z'^'ftY^+ps/o+z^fj/fl+z') * 
if- 1 ( p f (vr; 

Consider the indefinite integral 


(B. 15 ) 



A a 


ISrf(ltI l lAf+5 2 A 

l + ~ 


(B . 16 ) 
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This can be integrated by standard methods of integration; 


using the transformation 

A. A. 2 . -S' 2 


1 


___ \ izlX/i 

A A 

5, t -M 

and integrating, yields 


✓\ 2 . 

O 


(B . 17 ) 


£ 


1, = 




-ti it r 




■' hi mii L r+HYi'+icm/r) 

and, returning to the original variable, , 


18) 


£ 


I 


I . =7a-, ta / -4 


/\ 

f 

o 


? - 3 , 


3, 5/3. 1 

' 'I'" ‘l^l [d + f,?T+f+ 

Finally, by using Eq. (B.19), Eq. (b. 14) reduces to 

£ 

sf I 


I,m lit 


wi th 












f 7 ^ 


+ V { -#f ( f fc - ^ f ) 4- } 

1 III * 2.» fU J 


(B.19) 


(B.20) 


(B.21) 


67 


where 


2 = f f* + f 8 + ]' /2 

r t. > 

f%--riv+f; +?*:)* 

-[f'+r+fj ' 4 

( W>L L J W <"Vn 3 -* 


(B . 22 ) 


with 


^7*f ~ 


if Ja r ' "/r 

A A A A 

A A A £ 


(B. 23) 


= ?.» + M, 


A A 
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B.3 Source Distribution 

Finally, consider the integral of a source distribution 

of density a 5 /lySl over a trapezoidal planar element, given 
dX / 1 ' 1 
by 

£ _ fM _95 dZ 
I 5 = i/ R dX, |v ( S| 

(B . 2 4 ) 


, Y, P * p+dffYrY) 

__2^x 1 _ I dY f j_ ^ 

|aVa2,1 ^ i^r.-v 

where R is as shown in Eq. (B.8). Integrating Eq. (B.24) yields, 
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V- 


l uu v V ; + ]dl 


( i+* z tp 2 )' 


- f h M djtJ)* \[itAjArA l ) d f 1 


(B.25) 


A 
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A C C A A A 


where /> *7 . ^ , 5 f 1 t" are defined in Eqs. (B.15), 

u, / W it / >ov *j j 
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(B.13) and (B.9). Consider the indefinite integral 


£ 




A A A A 

with f — C * Integrating by parts yields (note 


A A 


e) “A ?. ? A + -|V) f- *(• 


Note that 

A; 

'l 


2 2. A A 

*1 i-e 


s-te ? ~ f-p ? 

. V 3 + ? 


Ai 

—1 _>f ) 

fV?‘ l f V 

I , f ^ 

■+ — + 


£ ' $vr ? 

. A A 

Combining Eqs . B.27 and B.28 yields (note that J ^ • 


1= * * i /-H - * 4 -fjftrf 


2* 


Note that 


ff d ? = f[ <i f M a? 2 V 3* 

*) 


I 0 ( * £>3, 


^1* -/FTP 

Combining Eqs. (B. 25) / (B. 29) and (b. 30) yields 


(B.26) 

that 

(B. 27) 


(B. 28) 



(B . 29 ) 


(B.30) 
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APPENDIX C 


GAUSSIAN QUADRATURE AROUND A SINGULAR POINT 


C,1 Introduction 

As mentioned in Subsection 4.3 f the coefficients 

are evaluated by approximating the surface element 2^ with 
its tangent plane at the centroid of the element. This value 
can eventually be corrected by adding the integral of difference 
between the original integrand and the (tangent-element) 
approximated integrand. This integration can be evaluated 
numerically by using standard Gaussian quadrature formulas. 
However, in the case of k = i (effect of the element on itself), 
the tangent plane contribution is equal to zero. Furthermore, 
the integrand becomes infinite when R^ = 0. Hence, a special 
integration scheme must be used. In this Appendix, an analysis 
of the type of singularity of the integrand of Eq. C. 1 (when 
k = i) is given (Subsection C.2). Then a transformation that 
eliminates the singularity is presented (Subsection C.3). 


C.2 Behavior of Doublets in the Neighborhood of Singularity 
For simplicity, the analysis of the behavior of the 
doublet in the neighborhood R = 0 is performed with a frame 
of reference such that the origin is at the centroid of the 
element and the Z-axis is directed as the normal n*. 
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Then, the equation of the element can be written as 

Z-F(X,r)^0 c .2 

with 2 = 0 for X = Y = 0 (Fig. C.l) and Eq. C.l for K = i, 
reduces to 

= -// dx < c ' 3 

where 

K= '§x,( x '' x) - §y;(1-Y) + (z,-z) c.4 

is the distance along the normal of the origin from the 
point X^, Y^ , (Fig. C.2). If R goes to zero, the distance 
h goes also to zero. 

k &R. Z C . 5 

where C = 1/ is the curvature of the cross section indi- 
cated in the figure ( is the radius of curvature). Thus, 

in a neighborhood of R = 0, Eq. C.3 reduces to 

C U = -jjr dX, dY, c.6 

It should be noted that C is the curvature of the cross section 
and thus C depends upon the angle of the cross section. Thus, 
in order to evaluate Eq. C.3, it is convenient to use polar 
coordinates since this eliminates the singularity 1/R as well 
as the sharp variations (in the plane X^, Y^) of the integrand, 
due to the dependence of C upon ^ . 
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<2.3 Integration Scheme 

As shown above, the integration of Eq. C.l (with k = i) 
in the neighborhood of the centroid of the element can be 
performed using standard quadrature techniques (Gaussian qua- 
drature in particular) in polar coordinates. However, with 
these variables, the definition of domain of integration is 
not given by coordinate lines. Hence, a more suitable tech- 
nique (fully correspondent to integration in polar coordinate) 
is described here. 

As shown in Eqs . 4.29 and 4.30, the boundary of the 
element in the plane are given by 

x c -Ax/ 2 s<x^y c +A//2 

y c - .< Y « Y c 

Note that the use of 

c * 

has the advantage of eliminating the square root singularity 
at the leading edge and the tip.* On the other hand, a 

singularity of the type R ^ in the plane X^, yields a 

singularity of the type , R ^ in the plane X, Y. Thus, the 

integral to be evaluated is of the type 


This singularity is due to the factor 


I VS) 


•'z-iskr' 


in 


■&S/JZ 
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where f (X, Y) is a finite but discontinuous function of X, Y 
(the discontinuity being due to the "cross-section-curvature 
effect") although continuous in polar coordinate. In order to 
analyze Eq. C.7 it is convenient to separate the contributions 
of + Dj and D 2 + (see Fig. C.3) as 

/ n C . 10 

C <N - C M + 

with 


c/. = // Hx,y)/(7\Y^ dxdr 

K* D x + 1> 3 

c; { • ff {(x.YiAr+Y'/'dxd? 

Using the transformation 

-I £ U£( 

Equation C. 11 reduces to 


C.ll 


C.12 


C . 13 


cjj = -ax^Y A' f(u,ir)luteludv- _ 
4 Jl v '('^ + ^WVfc+f’w/ 

Note that f(u, v) is a regular function of u 

~~2 ~2 - 1/2 

factor |u| compensates for the (X + Y ) ' 


ax AY, 


Tf^uv) 

-<A 


dud 1/ 

C . 14 


and v since the 
singularity and 
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the "cross-section-curvature effect" disappears in the u, v 
plane (which is similar to polar- coordinate plane) . Hence, 

Eq. C.14 can be evaluated by the Gaussian quadrature. Similar 
transformation can be used in Eq. €.12. 

This procedure was used to evaluate not only the effect 
Cj^ °f the element on itself, but also for effect of an ele- 
ment on the opposite element. Similar technique is used also 
for the evaluation of 

h = I, fa = £ £ % ^7 dl A C - 1: 
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X 


Fig. C.l Surface X. in neighborhood of R = 0 
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Fig. C.2 Curvature of cross section 
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